Attempt-time Monte Carlo: an alternative for simulation of stochastic jump processes 

with time-dependent transition rates 
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We present a new method for simulating Markovian jump processes with time-dependent tran- 
sitions rates, which avoids the transformation of random numbers by inverting time integrals over 
the rates. It relies on constructing a sequence of random time points from a homogeneous Poisson 
process, where the system under investigation attempts to change its state with certain probabili- 
ties. With respect to the underlying master equation the method corresponds to an exact formal 
solution in terms of a Dyson series. Different algorithms can be derived from the method and their 
power is demonstrated for a set of interacting two-level systems that are periodically driven by an 
external field. 
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I. INTRODUCTION 

Stochastic jump processes with time-dependent tran- 
sition rates are of general importance for many appli- 
cations in physics and chemistry, in particular for de- 
scribing the kinetics of chemical reactions [U-Q and the 
non-equilibrium dynamics of driven systems in statistical 
mechanics ■ With respect to applications in interdis- 
ciplinary fields they play an important role in connection 
with queuing theories. 

In general a system with N states is considered that at 
random time instants performs transitions from one state 
to another. In case of a Markovian jump dynamics the 
probability for the system to change its state in the time 
interval [t,t + At[ is independent of the history and given 
by Wij (t)At + o(At), where j and i ^ j are the initial and 
target state, respectively, and w%j(t) the corresponding 
transition rate at time t (wjj(t) = 0). This implies that, 
if the systems is in the state j at time to, it will stay in 
this state until a time t > to with probability <fij(t,to) — 
exp[— J* t * dr Wj 0t (r)], where w* ot (r) = J2i w ij( T ) i s tlie 
total escape rate from state j at time r. The probability 
to perform a transition to the target state i in the time 
interval [t,t + dt[ then is Wij(t)(j>j(t,to)dt, i. e. 



ipij(t,t ) = Wij(t)exp 



drwf^r) 



(1) 



is the probability density for the first transition to state i 
to occur at time t after the system was in state j at time 
to. Any algorithm that evolves the system according to 
Eq. JT]) generates stochastic trajectories with the correct 
path probabilities. 

The first algorithm of this kind was developed by Gille- 
spie 0] in generalization of the continuous-time Monte- 
Carlo algorithm introduced by Bortz et al. [1] for time- 
independent rates. We call it the reaction time algorithm 
(RTA) in the following. The RTA consists of drawing a 
random time t from the first transition time probabil- 



ity density ^ tot (Mo) = E^(Mo) = u;fVj(Mo) = 
—dt4>j{t,to) to any other state i ^ j, and a subsequent 
random selection of the target state i with probability 
Wij(t)/Wj 0t (t). In practice these two steps can be per- 
formed by generating two uncorrelated and uniformly 
distributed random numbers r%, T2 in the unit interval 
[0,1[ with some random number generator, where the 
first is used to specify the transition time t via 



W (t,to) 



dW ot (r) = -log(l-ri) (2) 



to 



and the second is used to select the target state i by 
requiring 



W k j(t) 



i-l 



< r 2 < 



E 



Wkj{t)_ 

(t) 



fc=1 < 



(3) 



Both steps, however, lead to some unpleasant problems 
in the practical realization. 

The first step according to Eq. ([2]) requires the calcu- 
lation of Wj(t,to) and the determination of its inverse 
Wj(.,t ) with respect to t in order to obtain the transi- 
tion time t = Wj(— log(l — r±),to)- While this is always 
possible, since w* ot > and accordingly Wj(t,to) is a 
monotonously increasing function of t, it can be CPU 
time consuming in case Wj(t,to) cannot be explicitly 
given in an analytical form and one needs to implement 
a root finding procedure. 

The second step according to Eq. ([3]) can be cumber- 
some in case there are many states (N large) and a sys- 
tematic grouping of the Wij (t) to only a few classes is not 
possible. This situation in particular applies to many- 
particle systems, where N typically grows exponentially 
with the number of particles, and the interactions (or 
a coupling to spatially inhomogeneous time-dependent 
external fields) can lead to a large number of different 
transitions rates. Moreover, even for systems with simple 
interactions (as, for example, Ising spin systems), where 



2 



a grouping is in principle possible, the subdivision of the 
unit interval underlying Eq. ([3]) cannot be strongly sim- 
plified for time-dependent rates. 

A way to circumvent Eq. is the use of the First Re- 
action Time Algorithm (FRTA) for time dependent rates 
H, or modifications of it 0]. In the FRTA one draws 
random first transition times tk from the probability den- 
sities ip k j(tk,t ) = w kj (t k )exp[- d,Tw k j(T)} for the 
individual transitions to each of the target states k and 
performs the transition i with the smallest ti — minj.{tj.} 
at time ti. This is statistically equivalent to the RTA, 
since for the given initial state j, the possible transitions 
to all target states are independent of each other. In 
short-range interacting systems, in particular, many of 
the random times tk can be kept for determining the 
next transition following i. In fact, all transitions from 
the new state i to target states k can be kept for which 
Wki{r) — Wkj(r) for r > t (see Ref. [l(| for details). 
However, the random times tk need to be drawn from 
"0fcj (tk, to) and this unfortunately involves the same prob- 
lems as discussed above in connection with Eq. @. 



II. ALGORITHMS 

We now present a new "attempt time algorithm" 
(ATA) that allows one to avoid the problems associated 
with the generation of the transition time in Eq. ([2]). 
Starting with the system in state j at time to as before, 
one first considers a large time interval T and determines 
a number /x*- ot satisfying 



> 



max 

t <T<t +T 



K ot (r)} . 



(4) 



In general this can by done easily, since Wj ot (T) is a 
known function. In particular for bounded transition 
rates it poses no difficulty, as, for example, in the case 
of Glauber rates or a periodic external driving, where 
T could be chosen as the time period. If an unlimited 
growth of u>j ot with time were present (an unphysical sit- 
uation for long times), T can be chosen self-consistently 
by requiring that the time t for the next transition to 
another state i ^ j (see below) must be smaller than 
to + T. 

Next an attempt time interval At\ is drawn from the 
exponential density Fj(Ati) = ^* ot exp(— ^* ot Ati) and 
the resulting attempt transition time t\ = to + At\ is 
rejected with probability pf 1 (t x ) = 1 - wf l (h) //z*. ot . If it 
is rejected, a further attempt time interval At2 is drawn 
from Fj(At2), corresponding to an attempt transition 
time ti = ti + A^2, and so on until an attempt time 
t < to + T is eventually accepted. Then a transition to 
a target state i is performed at time t with probability 
Wij (t) /Wj 0t (t) , using the target state selection of Eq. © . 

In order to show that this method yields the correct 
first transition probability density ipij(t, to) from Eq. (JTJ, 
let us first consider a sequence, where exactly n > at- 
tempts at some times t\ < . . . < t n are rejected and then 



the (n + l)th attempt leads to a transition to the target 
state i in the time interval [t, t + dt[. The corresponding 
probability density ip^j\t,to) is given by 



Vi n) (Mo)= / dt r , 

Jto Jto 



dtr 



dtt 



Wjj(t) 

w) ot (t) 



1 -pf (t)] Fj{t - t n ) J] pfit^Fjitm - tm-i) 



(5) 



n! 



t -I n 

dr^pf{r) 



Hf^t-ta)- / drwf\T) 



Summing over all possible n hence yields 



V>j(Mo) 



^2 ^if (*' <0 ) = Wi i (*) exp 



dTwf\T) 
(6) 

from Eq. (fTJ. 

It is clear that for avoiding the root finding of Eq. @ 
by use of the ATA, one has to pay the price for intro- 
ducing rejections. If the typical number of rejections can 
be kept small and an explicit analytical expression for t 
cannot be derived from Eq. ©, the ATA should become 
favorable in comparison to the RTA. Moreover, the ATA 
can be implemented in a software routine independent 
of the special form of the Wij (r) for applicants who are 
not interested to invest special thoughts on how to solve 
Eq. ©. 

One may object that the ATA still entails the problem 
connected with the cumbersome target state selection by 
Eq. However, as the RTA has the first reaction vari- 
ant FRTA, the ATA has a first attempt variant. In this 
first attempt time algorithm (FATA) one first determines, 
instead of /x* ot from Eq. upper bounds for the indi- 
vidual transitions to all target states k ^ j (fijj = 0), 



Mfcj > , max {w k j(T)} 

t <T<t + -l 



(7) 



Thereon random time intervals Atk are drawn from 
F k j(Atk) = fJ-kj exp(-/ijyAtfe), yielding corresponding 



attempt transition times t 



(i) 



Atk- The tran- 



sition to the target state k' with the minimal ijjy = 
mmfe-fij^} = t\ is attempted and rejected with probabil- 
ity ppkitff) = l-Wk'kitkh/Vk'k- If it is rejected, a fur- 

(2) (2) 

ther time interval At k , is drawn from Fk>j(At k , J ), yield- 
ing ty = t k X ) + At^\ while the other attempt transition 

times are kept, t^ — t^ for k =^ k' (it is not necessary 
to draw new time intervals for these target states due to 
the absence of memory in the Poisson process) . The tar- 

(2) (2) 

get state k" with the new minimal t k ,, — miiik{t k } = t% 
is then attempted and so on until eventually a transition 



3 



to a target state i is accepted at a time t < to + T. The 
determination of the minimal times can be done effec- 
tively by keeping an ordered stack of the attempt times. 
Furthermore, as in the FRTA, one can, after a successful 
transition to a target state i at time t, keep the (last up- 
dated) attempt times tk for all target states that are not 
affected by this transition (i. e. for which Wki(r) = Wkj(r) 
for t > t). Overall one can view the procedure implied 
by the FATA as that each state k has a next attempt 
time tk (with tj — oo if the system is in state j) and that 
the next attempt is made to the target state with the 
minimal tk- After each attempt, updates of some of the 
tk are made as described above in dependence of whether 
the attempt was rejected or accepted. 

In order to prove that the FATA gives the i()j(t,to) 
from Eq. |T]), we show that the probability den- 
sities Xij{t,t n ) = («)/«#* (t)](l - pfit))F 3 (t - 



t„) 

rej 



w ij (t)exp[-fif t (t - t„)] and %(t m , tm-i) 



tot 
3 



u>* ot ] exp[ 



-u tot (t 



p] el (t m )Fj(t m - t m -i) = [jit 

t m _i)] appearing in Eq. (J5J) are generated, if we set 
A** ot — J2k Mfcj (note that Eq. (Q| is automatically satis- 
fied by this choice) . These probability densities have the 
following meaning: Xij(P,t n )dt is the probability that, if 
the system is in state j at time t n , the next attempt to 
a target state occurs in the time interval [t,t + dt[, the 
attempt is accepted, and it changes the state from j to 
i'i f]j{t m ,t m -\)dt m is the probability that, after the at- 
tempt time t m , the next attempt occurs in [t m , t m + dt m [ 
with t m > t m _i and is rejected. 

In the FATA the probability Kij(t m ,t m —i)dt m that, 
when starting at time t m -i, the next attempt is occurring 
in [t m , t m + dt m [ to a target state I is given by 

d,Tfi kj exp(-fi kj T) 



Mj exp[-/i*. 0t (*m - *m-l)] 



(8) 



The product ensures that t m is the minimal time (the 
lower bound in the integral can be set equal to (t m — t m — i) 
for all k ^ I due to the absence of memory in the Poisson 
process). The probability that this attempted transition 
is rejected is p[J J (t m ) = 1 — wij(tm) / Hij and accordingly, 
by summing over all target states I, we obtain 



(t m ,t m -i) = 22p X if(tm)lJ>lj exp[-/i* ot )(i m - t m -i) 



tot „,,tot 



(tm)] exp[-^ ot (i m - t m -i)] 

(9) 



in agreement with the expression appearing in Eq. ([5]). 
Furthermore, when starting from time t n , the probability 
density X*i(^*n) referring to the joint probability that 
the next attempted transition occurs in [t,t + dt[ to state 



i and is accepted is given by 

Xij(t,t n ) = '^^-Kij(t,t n ) = w ij (t)&qp[-fif t (t - t n )] . 

(10) 

Hence one recovers the decomposition in Eq. ([5j with 

Before discussing an example, it is instructive to see 
how the ATA (and RTA) can be associated with a solu- 
tion of the underlying master equation 



^G(t,t') = -M(t)G(t,t') , 



G(t',t')=I (11) 



where G(t, t') is the matrix of transition probabilities 
Gij(t,t') for the system to be in state i at time t if it 
was in state j at time t' < t, and M(t) is the tran- 
sition rate matrix with elements Mj,j(t) — — my(i) for 
i ? j and Mji(t) = -Ei&Mi^t) = wf^t). Let 
us decompose M(t) as M(t) = D + A(t), where © = 
diag {^i ot , . . . , A^tv*}- If ^(*) were missing, the solu- 
tion of the master equation (TlTj) would be <&o(t,t') — 
diag {exp(-/4 ot (i - t'), exp(-(iffi(t - t')}. Hence, 
when introducing A(t,t') = G 1 (t, t')A(t)G (t, t') = 
Go(t', t)A(t)Go(t, t') in the "interaction picture", the so- 
lution of the master equation can be written as 



G(M') = Go(M') 



1+ / diiA(ii,i') 



dt 2 J eftiA(i 2 , *')■&(*!, + 



(12) 



Inserting 



_1 ]D> after each matrix A, one arrives at 



G(t,t') = G (t,t') + J <ft 1 G (t,*i)B(ti)fb(*i,*') 

+ / dt 2 ^ 2 ^iGo(t,i 2 )B(t 2 )Fo(^,ii)B(i 1 )Fo(ii,i') 

+ ... (13) 

where W (t,t') = BG (t,t') = diag-^ * exp[-fj,\ ot (t - 
t')], ffi exp[-fiffi(t - t')]}, and B(t) = A^)©" 1 has 
the matrix elements Bij(t) = —Wij(t) / '/i* ot for i ^ j and 



B n (t) = l 



fc (*)M 



Equation (1131) resembles the ATA: The transition prob- 
abilities Gij(t,t') are decomposed into paths with an ar- 
bitrary number n = 0,l,2,...of "Poisson points" , where 
transitions are attempted. The times between successive 
attempted transitions are exponentially distributed ac- 
cording to the matrix elements of Fo and the attempted 
transitions are accepted or rejected according to the 
probabilities encoded in the diagonal and non-diagonal 
elements of the B matrix, respectively. The Go enter- 
ing Eq. (|13p takes care that after the last attempt in a 
path with exactly n attempted transitions no further at- 
tempt occurs and the system remains in the target state 
i. The RTA can be associated with an analogous formal 
solution of the master equation if one replaces Gg(t, t') 
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by G* TA (M') = diag{«;J *(t)eocp[-/^dr«;i ot (T)],..., 
wffi(t)exp[- f* dTW^ir)]} and M(t) by B RTA (i) with 

elements £§ TA (t) = (1 - 5 l0 )w l0 {t)/wf i {t) (the diagonal 06 
elements are zero since the RTA is rejection- free). 

0.5 



III. EXAMPLE 

Let us now demonstrate the implementation of the 
FATA in an example. To this end we consider three 
mutually coupled two-level systems that are periodically 
driven. For an arbitrary given i, i = 1,2,3, the state 
\i, ± ) has the energy ±E(t). The occupancy of the state 
i,±) is specified by the occupation number rij = ±1. 
For example, if rii = —1, the 7-th two level system resides 
in the state | i, — ) and it possesses the energy —E(t). The 
coupling is described by the (positive) interaction param- 
eter V. The total energy of the three coupled two-level 
systems is given by the expression 



H(n, t) = V(n-in 2 + nxn 3 + n 2 n z ) + E(t) 



3 




(14) 



FIG. 1: Work distributions p(w) as obtained from the FATA 
for AE = V = 5 and uj = 0.1 (squares, green color), 1 (circles, 
blue color), and cj = 10 (stars, red color). 



where n = (m, 712,713) specifies the microstate of the 
compound system. The periodic driving is considered 
to change energies of the individual two-level systems as 



E(t) = — sin(wt) , 



(15) 



where AE > is the amplitude of modulation and uj 
its frequency. Due to contact of the compound system 
with a heat reservoir at temperature T, transitions be- 
tween its microstates occur. Assume that in the initial 
state m one and only one occupation number differs from 
the corresponding occupation number in the final state 
n. Then instantaneous value of the detailed-balanced 
Glauber jump rates connecting these two states reads 



w(m — > n, t) 



l + exp{P[H(n,t)-H(m,t)}} 



(16) 



The other pairs of microstates are not connected, that is, 
the transition rates between them vanish. In the above 
expression, v designates an attempt frequency, and /3 is 
the inverse temperature. In the following we will use k^T 
as our energy unit and v~ l as our time unit. 

In current research of non-equilibrium systems, in par- 
ticular of processes in small molecular systems, the in- 
vestigation of distributions of microscopic work receives 
much attention. Among others, this is largely motivated 
by questions concerning the optimization of processes, 
and by the connection of the work distributions to fluc- 
tuation theorems. These theorems allow one to obtain 
equilibrium thermodynamic quantities from the study of 
non-equilibrium processes and they are useful for getting 
a deeper insight into the manifestation of the second law 
of thermodynamics. At the same time, the analytical ex- 
pressions for the work distribution are rarely attainable 



(one exception is reported in |l2|). It is therefore inter- 
esting to see how the FATA can be employed for studies 
in this research field. To be specific, we focus on the sta- 
tionary state and calculate work distributions within one 
period of the external driving. For these distributions we 
check the detailed fluctuation theorem of Crooks [IJ, as 
generalized by Hatano and Sasa [ll| to steady states (for 
a nice summary of different forms of detailed and integral 
fluctuation theorems, see @). 

In our model, due to the possibility of thermally acti- 
vated transitions between the eight microstates, the state 
vector n must be understood as a stochastic process. We 
designate it as n(t), and let n tr (i) denotes its arbitrary 
fixed realization. The instantaneous energy of the com- 
pound system along this realization is then H(n tr (t), t). 
The work done on the system during the mth period 
[hit, (771 + 1)t], t = 2t:/uj, if the system evolves along the 
realization in question, is given by 



(m+l)r p, 
rnr ^ 



jAE 



E 



(m+l)r 



dtnf(t)cos(ut) . (17) 



In the stationary limit m — > 00 (m S> 1) we can drop 
the index m. According to the detailed fluctuation theo- 
rem, the work distribution p(w) should, in our case (time- 
symmetric situation with respect to the initial microstate 
distribution for starting forward and backward paths), 
obey the relation p(w) exp(-w) = p{— w). 

Figure [T] shows the results for p(w) obtained from the 
FATA for AE = V = 5, and three different frequencies 
uj = 0.1, 1, and 10. First, we let the system evolve during 



5 




0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 

p(-ui) 



FIG. 2: Check of the detailed fluctuation theorem 
p(w) exp(-w) = p(—w) for the work distributions shown in 
Fig. \T\ The same symbols/colors are used for the three dif- 
ferent frequencies as in Fig. [1] 



the N ini = 1 (w = 0.1), N- mi = 3 (w = 1), iV lnl = 9 (w = 
10) periods to reach the stationary state. Subsequently, 
the work values w tr according to Eq. (|T7|) were sampled 
over N = 10 4 (oj = 0.1), N = 10 5 (oj = 1), and N = 10 5 
(to = 10) periods. 

With decreasing cj, the maxima of the work distribu- 
tions in Fig. Q] shift toward w = 0, and ^-singularities, 
marked by the vertical lines, receive less weight. These 
5— singularities are associated with stochastic trajectories 
of the system, where no transitions occur within a period 
of the driving. For oj = 0.1, p(w) is already close to the 
Gaussian fluctuation regime. 

In Fig. [5] we show that the work distributions from 
Fig.[T]indeed fulfill the detailed fluctuation theorem. This 
demonstrates that the FATA successfully generates sys- 



tem trajectories with the correct statistics of the stochas- 
tic process. 



IV. SUMMARY 

In summary, we have presented new simulation al- 
gorithms for Markovian jump processes with time- 
dependent transition rates, which avoid the often cum- 
bersome or unhandy calculation of inverse functions. The 
ATA and FATA rely on the construction of a series of 
Poisson points, where transitions are attempted and re- 
jected with certain probabilities. As a consequence, both 
algorithms are easy to implement, and their efficiency will 
be good as long as the number of rejections can be kept 
small. For complex interacting systems, the FATA has 
the same merits as the FRTA with respect to the FRA. 
Both the ATA and FATA generate exact realizations of 
the stochastic process. Their connection to perturbativc 
solutions of the underlying master equation may allow 
one to include in future work also non-Markovian fea- 
tures of a stochastic dynamics by letting the rejection 
probabilities to depend on the history |13| . Compared to 
the RTA and FRTA, the new algorithms should in par- 
ticular be favorable, when considering periodically driven 
systems with interactions. Such systems are of much cur- 
rent interest in the study of non-equilibrium stationary 
states and we thus hope that our findings will help to 
investigate them more conveniently and efficiently. 
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